#Project: TV polarization
#Author: Josh Ryan and Jeff Lyons
#Working Date: Winter 23
#Task: Use Imai, Kim, Wang PanelMatch analysis--analysis using ps match--for legislator level, Figure 9
#Related files: leg_level_forstata_v6_102524.dta and chamber_level_forstata_v2_112623.dta is loaded
#Updated for PanelMatch 3.1.1

library(foreign)
library(effects)
library(sciplot)
library(sandwich)
library(lmtest)   
library(BMS)
library(rJava)
library(lattice)
library(plyr)
library(lme4)
library(gmodels)
library(MASS)
library(stringr)
library("XLConnect")
library(reshape2)
library(xlsx)
library(dplyr)
library(doBy)
library(Hmisc)
library(stargazer)
library(xtable)
library(dplyr)
library(xtable)
library(foreign)
library(dotwhisker)#https://cran.r-project.org/web/packages/dotwhisker/vignettes/dotwhisker-vignette.html
library(berryFunctions) #for the insert row function
library(ggplot2)
library(gridExtra)
library(wnominate)
library(tidyr)
library(wnominate)
library(oc)
library(haven)
library(PanelMatch)

dat<-read_dta("leg_level_forstata_v6_102524.dta")

dat.cham<-read_dta("chamber_level_forstata_v4_102524.dta")

####Figure I6

dat.cham<-dat.cham[c("state", "chamber", "state_cham")]
dat.cham$state_cham<-as.numeric(dat.cham$state_cham)
dat.cham$chamber<-tolower(dat.cham$chamber)

dat.cham<-dat.cham %>%
  distinct(state_cham, state, chamber, .keep_all = TRUE)

dat$state_cham<-NULL
dat<-merge(dat, dat.cham, by=c("state", "chamber"), all.x = TRUE)

dat<-subset(dat, state_cham==71 | state_cham==72 | state_cham==43 | state_cham==44 | 
              state_cham==57 | state_cham==58 | state_cham==49 |state_cham==50 | 
              state_cham==27 | state_cham==28 | state_cham==5 | state_cham==6 | 
              state_cham==47 | state_cham==48 | state_cham==73 | state_cham==74 | 
              state_cham==39 | state_cham==40 | state_cham==61 | state_cham==62 | 
              state_cham==7 | state_cham==8 |state_cham==51 | state_cham==52 | 
              state_cham==96 | state_cham==54 | state_cham==79 | state_cham==80 |
              state_cham==23 | state_cham==24 | state_cham==69 | state_cham==70 | 
              state_cham==75 | state_cham==76 | state_cham==35 | state_cham==91 | 
              state_cham==92 | state_cham==87 |state_cham==88)

#Have to keep making as dataframe because of the way Panelmatch works
dat<-as.data.frame(dat)

#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.np<-dat[c("np_score", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.np<-as.data.frame(dat.np)
dat.np$distid1<-as.integer(dat.np$distid1)
dat.np$year<-as.integer(dat.np$year)

dat.np0<-subset(dat.np, party_id==0)
dat.np1<-subset(dat.np, party_id==1)

#check balance improvement first run matching results, then look at balance improvement
######np score Dems

# Create PanelData object for Democrats
panel_data_dem <- PanelData(
  panel.data = dat.np0,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "np_score"
)

results.mah.ideo.dem <- PanelMatch(
  panel.data = panel_data_dem,
  lag = 4,
  refinement.method = "mahalanobis",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.ideo.dem<-results.mah.ideo.dem$att
PE.results.ideo.dem <- PanelEstimate(sets = results.mah.ideo.dem, panel.data = panel_data_dem)
summary(PE.results.ideo.dem)
PE.results.ideo.dem<-summary(PE.results.ideo.dem) #save as dataframe
names(PE.results.ideo.dem)
coef<-PE.results.ideo.dem[ , 1] #get time coefficients
se<-PE.results.ideo.dem[ , 2] #get standard errors
lower<-PE.results.ideo.dem[ , 3] #lower confidence bound
upper<-PE.results.ideo.dem[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

results.ideo.dem<-as.data.frame(cbind(coef,se,lower,upper,time))
results.ideo.dem

plot.ideo.dem<-ggplot(results.ideo.dem, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-.2, .2, .04)) +
  coord_cartesian(ylim=c(-.2, .2)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Democratic Ideology")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.ideo.dem)

######################GOP
# Create PanelData object for Republicans
panel_data_gop <- PanelData(
  panel.data = dat.np1,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "np_score"
)

#covariate balance propensity score matching
results.cbps.ideo.gop <- PanelMatch(
  panel.data = panel_data_gop,
  lag = 4,
  refinement.method = "CBPS.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.ideo.gop<-results.cbps.ideo.gop$att
PE.results.ideo.gop <- PanelEstimate(sets = results.cbps.ideo.gop, panel.data = panel_data_gop)
summary(PE.results.ideo.gop)
PE.results.ideo.gop<-summary(PE.results.ideo.gop) #save as dataframe
names(PE.results.ideo.gop)
coef<-PE.results.ideo.gop[ , 1] #get time coefficients
se<-PE.results.ideo.gop[ , 2] #get standard errors
lower<-PE.results.ideo.gop[ , 3] #lower confidence bound
upper<-PE.results.ideo.gop[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

results.ideo.gop<-as.data.frame(cbind(coef,se,lower,upper,time))
results.ideo.gop

plot.ideo.gop<-ggplot(results.ideo.gop, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-.2, .2, .04)) +
  coord_cartesian(ylim=c(-.2, .2)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Republican Ideology")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.ideo.gop)

###############SLES

#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.sles<-dat[c("sles", "party_id", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

#have to constantly be making these as integers, super annoying
dat.sles<-as.data.frame(dat.sles)
dat.sles$distid1<-as.integer(dat.sles$distid1)
dat.sles$year<-as.integer(dat.sles$year)

# Create PanelData object for SLES
panel_data_sles <- PanelData(
  panel.data = dat.sles,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "sles"
)

results.ps.sles <- PanelMatch(
  panel.data = panel_data_sles,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4))+ I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + 
    I(lag(mds1, 1:4)) + I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.sles<-results.ps.sles$att
PE.results.sles <- PanelEstimate(sets = results.ps.sles, panel.data = panel_data_sles)
summary(PE.results.sles)
PE.results.sles<-summary(PE.results.sles) #save as dataframe
names(PE.results.sles)
coef<-PE.results.sles[ , 1] #get time coefficients
se<-PE.results.sles[ , 2] #get standard errors
lower<-PE.results.sles[ , 3] #lower confidence bound
upper<-PE.results.sles[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

PE.results.sles<-as.data.frame(cbind(coef,se,lower,upper,time))
PE.results.sles

plot.sles<-ggplot(PE.results.sles, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-1, 1, .4)) +
  coord_cartesian(ylim=c(-1, 1)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Legislative Effectiveness")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.sles)

###############Party Loyalty--Too many missing to work

# Uncomment and use if you want to include party loyalty analysis
# #keep only relevant variables, doesn't matter but panelmatch throws an error
# dat.loyalty<-dat[c("party_loyalty", "party_id", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]
# 
# #have to constantly be making these as integers, super annoying
# dat.loyalty<-as.data.frame(dat.loyalty)
# dat.loyalty$distid1<-as.integer(dat.loyalty$distid1)
# dat.loyalty$year<-as.integer(dat.loyalty$year)
# 
# # Create PanelData object for Party Loyalty
# panel_data_loyalty <- PanelData(
#   panel.data = dat.loyalty,
#   unit.id = "distid1",
#   time.id = "year",
#   treatment = "treatment",
#   outcome = "party_loyalty"
# )
# 
# results.psmatch.loyalty <- PanelMatch(
#   panel.data = panel_data_loyalty,
#   lag = 4,
#   refinement.method = "ps.match",
#   match.missing = TRUE,
#   covs.formula = ~ I(lag(majority, 1:4))+ I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + 
#     I(lag(mds1, 1:4)) + I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
#     I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#     I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#   size.match = 10,
#   qoi = "att",
#   lead = 0:4,
#   forbid.treatment.reversal = TRUE
# )
# 
# mset.loyalty<-results.psmatch.loyalty$att
# PE.results.loyalty <- PanelEstimate(sets = results.psmatch.loyalty, panel.data = panel_data_loyalty)
# summary(PE.results.loyalty)
# PE.results.loyalty<-summary(PE.results.loyalty) #save as dataframe
# names(PE.results.loyalty)
# coef<-PE.results.loyalty[ , 1] #get time coefficients
# se<-PE.results.loyalty[ , 2] #get standard errors
# lower<-PE.results.loyalty[ , 3] #lower confidence bound
# upper<-PE.results.loyalty[ , 4] #upper confidence bound
# time<-c(0,1,2,3,4)
# 
# PE.results.loyalty<-as.data.frame(cbind(coef,se,lower,upper,time))
# PE.results.loyalty
# 
# plot.loyalty<-ggplot(PE.results.loyalty, aes(time,coef)) +
#   geom_point() +
#   geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
#   theme_classic() +
#   scale_y_continuous(breaks=seq(-.2, .2, .04)) +
#   coord_cartesian(ylim=c(-.2, .2)) +
#   scale_x_continuous(breaks=seq(0, 4, 1)) +
#   theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank(), 
#         panel.background = element_blank(),
#         axis.line = element_line(colour = "black"))+
#   theme(axis.line.x = element_line(color="black", size = .5),
#         axis.line.y = element_line(color="black", size = .5))+
#   geom_hline(yintercept=0, color="gray")+
#   xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
#   ggtitle("Party Loyalty")+theme(plot.title = element_text(hjust = 0.5))
# 
# plot(plot.loyalty)

all.plots<-grid.arrange(plot.ideo.dem, plot.ideo.gop, plot.sles, ncol=2)

plot(all.plots)


